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S^ ■ Abstract 
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ti . The recent experimental discoveries about excitation energy transfer (EET) in hght harvesting 

-*— > \ 

C^ , antenna (LHA) attract a lot of interest. As an open non-equilibrium quantum system, the EET 

c/3 ■ 

) ' ' demands more rigorous theoretical framework to understand the interaction between system and 

a : 

^ ' environment and therein the evolution of reduced density matrix. A phonon is often used to 

'^ , model the fluctuating environment and convolutes the reduced quantum system temporarily. In 

O I this paper, we propose a novel way to construct complex-valued Gaussian processes to describe 

thermal quantum phonon bath exactly by converting the convolution of influence functional into 

(n: 

K* . the time correlation of complex Gaussian random field. Based on the construction, we propose a 

m ■ 

'^ \ rigorous and efficient computational method, the covariance decomposition (CD) and conditional 

^^ ' propagation scheme, to simulate the temporarily entangled reduced system. 

^: 

^^ ■ The new method allows us to study the non-Markovian effect without perturbation under the 

m ; 

'""' , influence of different spectral densities of the linear system-phonon coupling coefficients. Its ap- 

"^ ' plication in the study of EET in the Fenna-Matthews-Olson (FMO) model Hamiltonian under 

H ' four different spectral densities is discussed. Since the scaling of our algorithm is linear due to its 

Monte Carlo nature, the future application of the method for large LHA systems is attractive. In 

addition, this method can be used to study the effect of correlated initial condition on the reduced 

dynamics in the future. 



I. INTRODUCTION 

The study of irreversible open quantum dissipative processes is important in almost every 
field of condensed-matter physics and chemistry, such as reaction rate theory, ultrafast 
phenomena, tunneling at defects in solids, and quantum optics, etc.—''—. The 2D spectroscopic 
experiments^ in the light harvesting antennas (LHA) reveal the existence of the long-last 
coherence in EET. As a result, the real-time dynamics of EEA in the photosynthetic light 
harvesting environment attracts a lot of theoretical interest and debates^. The almost perfect 
EET efficiency in the disordered LHA protein environment could be related to the preserved 
quantum coherence. As an open quantum system, the unique decoherence and relaxation 
due to the LHA protein environment should play an important role in EET-. 

Dynamic relaxation and decoherence due to environmental fluctuations contains the crit- 
ical information of reduced system dynamics and the interaction between system and bath. 
For classical systems with linear dissipation, Langevin equations (or Ito stochastic differential 
equation(SDE))^i^ and its extension provide a simple theoretical (and numerical) framework 
to describe the interaction between a system and a complex thermal reservoir in terms of 
stochastic forces and memory friction. For open quantum systems, the corresponding ac- 
count of quantum noise is still an open question. Developing a new and rigorous numerical 
methodology to simulate open quantum dynamics will provide a novel understanding of 
EET in biological processes. 

System bath Hamiltonian with bilinear coupling is the common model to study dissipative 
dynamics. Particularly, the exciton-phonon coupling Hamiltonian has been used to study 
exciton transport. In the Liouville space, the collective motion of environmental phonon 
modes is reduced to influence functional, which convolutes the reduced system dynamics. 
The convolution entangles reduced quantum dynamics. As a result, calculation complexity 
grows exponentially over the time. The traditional way to avoid the computational issue is 
using the different truncation schemes based on cumulant expansion or project operator-""— 
in the weak coupling limit. Therefore, the model based on second-order master equations can 
not accurately describe the interaction of system and bath and introduce the inconsistence 
about reduced system dynamics. Under some special cases, the integral equation based on 
influence functional can be reduced to a non-perturbative hierarchical equation of motion 
(HEOM)i^ii^ and similar ideas have been explored in different contexts^. But the HEOM 



method is very much hmited to the Drude-Lorentzian spectral density {i.e. exponential 
kernel.) Therefore, HEOM can not be used for arbitrary spectral densities, requires ad 
hoc cutoff of the infinite iterative equations and the scaling of the method is nonlinear and 
unclear. In addition, approaches based on semiclassical path integraU^"— , hybrid Ehrenfest 
and NIBA method^, iterative tensor product methoci^'^S and quantum Brownian motion 
process2ii2^ have been proposed recently and potentially can be used for the general spectral 
density cases. 

The convolution due to influence functional^^ makes the computation of reduced system 
dynamics extremely complicated. Due to the similarity between generalized characteristic 
function and influence functional, we can map the convolution kernel of influence functional 
to the covariance matrix of a Gaussian process. 

Convolution i — )■ Correlation. (1) 

So, we can linearize the computational effort with the covariance decomposition (CD) 
method to sample the quantum fluctuation of a phonon environment as a complex-valued 
Gaussian random field. In other words, constructing Gaussian random field with the CD 
method can deconvolute the reduced system dynamics. This method will give us a general 
computational tool to study the effect of phonon environment on reduced system dynamics 
in a rigorous and complete way. The sampling strategy for complex-valued Gaussian random 
process of arbitrary temperature will be more complicated than the real Gaussian random 
process at the high temperature limit. We will show how to construct the complex- valued 
Gaussian random process in this paper. However, I will present the computational results 
for the high temperature limit, i.e., ignoring the imaginary part of the kernel and quantum 
detailed balance^^. The sampling strategy for the complex-valued Gaussian random process 
based on the complex unitary transformation will be the future work. 

The coupling part of system bath Hamiltonian is critical in determining the interaction 
between system and environment. Since the coupling is bilinear in the exciton-phonon 
coupling Hamiltonian, the spectral density of the coupling strength between the system and 
bath is the major factor in relaxation and decoherence processes^^. The CD method will 
be a computationally efficient tool to study the effect of the arbitrary distribution (spectral 
density) of the linear coupling coefficients on the evolution of reduced density matrix. Two 
major challenges from the interaction in the current research work are : 1. the quantum 



memory effect of plionon environment (Non-Markovianity); and 2. tlie correlated initial 
conditions (inseparability of the system and bath). The CD method will allow us to address 
the two challenges in the same framework. Currently, the CD method is base on the influence 
functional formalism so that the original assumptions of influence functional will be kept^^. 
Some authorsi^i2i"— have used different approaches to construct to complex- valued Gaussian 
processes. But the approach we offer here is based on multivariate complex-valued normal 
distribution function and can take advantage of the established numerical Monte Carlo 
methods based on normal distribution functions. Our approach is indifferent to the kernel of 
influence functional and the choice of spectral density. For example, you can only generate 
a Markov chain for a Gaussian process with a exponential kernel (Gauss-Markov model) 
corresponding to the high temperature limit of Drude spectral density. For the general 
spectral density, it is impossible to generate a Markov chain to sample the corresponding 
Gaussian process. So our method is general and efficient. 

In the discussion section, we like to find out how the shape, such as slope, tail , and center 
of spectral densities can change the relaxation and coherence of reduced system dynamics. 
It is particularly interesting to examine the optical phonon band that have different spectral 
density from the acoustic one^^. We will look at the geometric impact of the spectral density 
in this paper. 

The paper is organized into five sections: 1. in Sec. Ull we briefiy review coherent state 
path integral and infiuence functional formalism; 2. in Sec. IIIIl we introduce the generalized 
characteristic function and random evolution operator. With them, we drive the stochastic 
integral and differential equations. In addition, we discuss how to derive the Gauss-Markov 
model in our framework at the high temperature limit; 3. in Sec. llVt we introduce the CD 
method and conditional propagation scheme; 4. in Sec. |Vl we present the benchmarking of 
the conditional propagation scheme according to the Gauss-Markov model. We also show 
and discuss the results for FMO under the influence of different spectral densities. 

To go beyond the influence functional (the bilinear coupling of system and bath) for the 
reduced system dynamics, the path integral is the last resort, which will allow us to study 
the nonlinear coupling of the system and bath motions. 



II. INFLUENCE FUNCTIONAL 

In this section, we give a brief description of influence functional based on the following 
exciton-phonon coupling Hamiltoniani^i^, 

H = Hs{a\ a) + Hi{a\ a, X) + Hb, (2) 

where Hs{a), a) = Ylk ^k ^I Ofc, where a^ and a are the set of a\ and Ofc, Hi = V{a\ a) x X, 
X = ^jCjXj and Hb = Xlj ( 2m~ '^ |"^j'^|^i)- ^^^ interaction between system and bath 
are bilinear which decides that there is only one quanta of energy can be move in or out of 
the system every time. 

The evolution of the isolated system density matrix ps{tf) under the Hamiltonian of 
system Hs{a\ a) in the path integral representation can be expressed, 

psiz},z'j,tj) = \<ilitj)){<il'itf)\= (3) 

N~^dz*fdzf j j N~^dz'j*dz'j^\zf)K{z*j, z'j^,tf,ti){z'j\, 



where z and z' are complex c-number for bosons and its conjugate, and the kernel defined 

as, 

K{z},z'f,tf,U) = j Df[Q{T)]j Df[Q'{T)] (4) 

xexp[(z/n)55(Q,t/,t.)]exp[-(z/;i)5^(Q',t^,t,)], 

where the prime sign ' is the indicator of the coordinates associated with the 
bra state (vl/'(t^)|, /D/[Q(r)] = \imN^^Y.ti ! N-'dz*{n)dz{n), I DflQ'ir)] = 
lim7v->>oo X]j=i j N~^dz'*{Ti)dz'{Ti), N is the normalization factor of coherent states. 
Ss{Q,tf,ti) and Sg{Q',tf,ti) are actions defined as 

SsiQ, tf, U)= r dr {ih z*{t)z{t) - Hs{z*{r), z{r))) , (5) 

Jti 

and 

SUQ',tf,U) = r dr{-zh z'*{T)z'{r) + Hs{z'*{t),z'{t))) , (6) 

Jti 

Q{t) is the short notation for the pair of (z(r), 2;*(r)), Q'{t) for {z' {r) , z'* {t)) . Once the 
bath is coupled to the system, assuming that the initial total density matrix is separable 
Ptoi(O) = ps'(O) *p^(0), where pi is the equilibrium density matrix of the bath, -priac'i-BH )) ' 



the evolution of the reduced system psitf) = Tr(p(t)) can be expressed in the similar 
expression in Eq. Hlwith the new K{zf*, z'jr,tf,ti), 

K{z},z'f,tf,U) = j Df[Q{r)]j Df[Q\r)] (7) 

xexp[(V/i)5,(Q,t/,t,)] X F{Q,Q'-tf,U) x exp[-{i/h)S:{Q',tf,ti)], 

in which F{Q,Q' ,tf,ti), influence functional, is deflned as, 

F{Q, Q- tf, U) = exp{- [' dr [ da (l^(Q(r)) - V{Q'{t)) (8) 

[7(r - a)ViQia)) - ^ (r - a)V (Q' (a))]} , 

where 7 (t) = Li — iL2, Li(t) = J^ dujJ{uj) coth.{(3huj/2) cos{ujt), 1^2 (t) = J^ duJ{uj)sm{ujt), 
and J{co) = 2 Ylj m^ui- ^^^ ~ '^i) ^^ ^^^ spectral density, which describes the distribution of 
the coupling strength coefficients between the system and different Harmonic modes. For 
the simplification, we assume h = 1 from now on. 

In the next section, we will show how to map the convolution of infiuence functional to 
the correlation of Gaussian random process. 

III. GENERAL CHARACTERISTIC FUNCTION AND RANDOM EVOLUTION 
OPERATOR 

Our construction of Gaussian random process mathematically is based on the general 
characteristics function of classical discrete Gaussian process (DGP). The details of the 
construction of DGP and general characteristics function (GCF) are briefiy reviewed in 
App. [XI By extending this construction, we can obtain the mapping defined in Eq. [1], i.e., 
reproducing the convolution of a reduced system quantum dynamics with a Gaussian random 
process and associated random evolution operator. The goal of this construction is to replace 
the convoluted path integral defined in Eq. [3] with the path integral conditional on the 
environment fiuct nation. 

A. Random Evolution Operator 

In order to map infiuence functional, we need to extend the real-valued covariance ma- 
trix in Eq. IA3I to a complex- valued covariance matrix with kernel, 7(r, a) in Eq. [HI As a 
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result, the complex-valued GCF will be equivalent to influence functional defined in Eq. |8l 
To accommodate the structure of double path integral in Eq. [3|, we propose the following 
complex-valued Gaussian stochastic process ^(t) = [^(t),^'{t)]. With the stochastic process, 
the convolution due to the influence functional will be decomposed to a random evolution 
operator. 



F{Q,Q'-tf,U) = {exp 



—I 



dr{V{Q{T))i{r)-V{Q\T))ar)) 



(9) 

m 



where ()iu\ is the expectation average over the trajectories ^(t), which are independent of 
Q{t) and Q'{t), therefore V{Q{t) and V{Q'{t). The detailed derivation can be found in 
App. IHl Since the Gaussian random process $,{t) is independent of the system operators, the 
reduced system density, ps in terms of path integral in Eq. [31 can be re-written as. 



Psitf) = j Df[Q{r)]j Df[Q\r)] (10) 

exp[(2/n)5,(Q(r),t/,t,)] x F{Q,Q';tf,U) x exp[-{i/h)S:{Q'{T),tf,U)] 



exp [{i/h){Ss{Q{T),tj,U) + V{Q{r)) * i{r)] 
exp [-{t/h){Sm{r).tf,U) + V{Q\t)) * ^'(r))] > (11) 



r exp S.-tf' dT[Hs + V{a^, a) ^(r)] | ^^(0) 



u 
rexpUj^'dr[Hs + Via\a)^'{r)]\y , (12) 

where T is time ordering operator. Eq. [T2] gives the linear random evolution operators for 
the forward and backward propagation. By sampling the trajectories of ^(t) and applying 
the above equation of motion, we can calculate the evolution of the reduce system density 
matrix, psit) = {p(t\i))^,^y With Eq. |9l we convert the influence functional convolution to 
an expectation average of the reduced system dynamics conditional on the Gaussian random 
fleld trajectories. The equivalent differential form of Eq. [12] can be expressed as, 

^P^ = -I Ls psm - I [V{a\ a) i{t) psiAl) - psiAl) V{a\ a) ^(t)], (13) 

where Ls = [Hs, ■]■ At the high temperature limit, the two distinct processes, ^(t) and ^'(t) 
will collapsed to one C,{t) and the term V{a'', a) C,{t) ps{t) — ps{t) ^(a^ a) C,'{t) will become 



a commutator bracket, [V^(a''',a) ^(t), p(t)]. Then we can recover the stochastic Liouville 
equation which is extensively used in the study of the exciton transport^, 

^ = -< m ,M), (14) 

where L{t) = [Hs + V{a\ a) ^{t), ■]. The details of the reduction of quantum noise to classical 
noise are discussed in App O 

The structure of influence functional determines that covariance matrix is non-Hermitian, 
so that it is different from the one proposed by Miller and Coworkers^^ for the signal pro- 
cessing and others. The difference is clearly reflected in Eq. IB3I Simply speaking, C,(t) and 
^'{t) are not conjugate to each other, which is the nature of the open quantum dynamics 
embed in the influence functional. 

B. Multichromophore Frenkel- Exciton System 

For the multichromophore Frenkel-Exciton system, we need to consider the path inter- 
ference. As a result, extra steps should be taken to replicate the right influence functional 
kernel since besides the time correlation, the correlation between different sites can exist 
due to path interference. Here, we take a dimer system as an example to explain the extra 
steps. The exciton-phonon coupling in a dimer is defined as, 

H = Hs + VxX + I xHb, (15) 

where 

Hs=\ ' ' I, (16) 




^h''°I, (17) 

V2 

X is the bath operator, and / is identity matrix. Therefore, in order to replicate the 
convolution kernel, we need two independent Gaussian random processes, ^i(t) and ^2(^) 
with the same kernel. The details of the construction of the Gaussian processes for the 
dimer Hamiltonian are presented in App. |Dl The independence means there is no spatial 
correlation between ^i(t) and ^2(^) o^ ^^e two Gaussian processes have two independent 
covariance matrix as defined in Eq. IB3I At the high temperature limit, if the phonon band 
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is only coupled to the site energy, the Hamiltonian of dimer in Eq. [15] can be reduced to the 
stochastic Hamiltonian, 

.„J--fW ^ V (18) 

If V has the off-diagonal matrix elements, Vtj, then spatial correlation will appear between 
different sites due to the interference between different paths as discussed in the paper of 
the enhanced coherence^. 

C. High Temperature Limit and Gauss-Markov Model 

Gauss-Markov (Ornstein-Uhlenbeck) model has been used for the study of exciton trans- 
port process with memory. The environmental fluctuation is modeled with the Ornstein- 
Uhlenbeck (OU) process which is a real Gaussian process with exponential kernel. Although 
Gauss-Markov model is a phenomenological model, we can derive the Gauss-Markov from 
Eq. [13] with a proper assumption of the spectral density. 

If the spectral density is a Drude-form Lorentzian function, J{co) = -^ 2T^2 , then Li = 
A^^exp(— 7((T — r)), and L2 = A^exp(— 7(0" — r)) as discussed in the literature^i^. It 
can be clearly shown that at the high temperature limit {i.e., (3 — > 0), ^ — t- and Li — ?■ 
0. However, L2 has the exponential form. As a result, the Gauss-Markov model can be 
considered as the high temperature limit of the complex- valued Gaussian process with the 
Drude-form Lorentzian spectral density. 

The Gauss-Markov model^^i^ has the stochastic Hamiltonian, H = Hg + C,{t), where ^(t) 
is the OU process with the kernel 7(r — a) = A^exp(cr — r). As a result, {C,{t)^{(t)) = 
{^'i'r)C{(T)) = (f (r)^(a)) = (^(r)f (a)) = 7(r - a) in Eq[Ba In other words, the Gauss- 
Markov model is a specific case of Eq. [T4]when the the Gaussian process kernel is exponential. 
The details of the reduction of complex-valued Gaussian processes to real-valued ones are 
presented in App. [C] 

IV. METHODOLOGY 

In this section, we propose the simulation method based on the derivations of Gaussian 
processes. Firstly, we will discuss the covariance decomposition (CD) method and Monte 



Carlo sampling strategy to generate the trajectories of environment fluctuations. Secondly, 
we present the conditional propagation scheme to simulate the reduced system quantum 
dynamics conditional on the trajectories generated from the CD method. 

A. Covariance Decomposition and Monte Carlo Sampling of Environment Fluc- 
tuation 

The Gaussian random process of our construction is independent of the system degree of 
freedom. Therefore we can sample them before we calculate the reduced system propagation 
conditional on the Gaussian random process. With the complex-valued PDF in Eq. IB5[ 
potentially we can sample the discrete complex- valued Gaussian process on the discrete time 
lattice. At the high temperature, we show in App. [Cl the complex-valued PDF becomes 
a real-valued one. For the real-valued PDF, the Cholesky decomposition (CD)^ is used 
to sample the discrete Gaussian trajectories by transforming correlated random variables 
to uncorrelated ones, which decomposes the covariance matrix in Eq. ICll into a product of 
lower and upper triangle matrix, A = LL?" where A = 7^^, the inverse of covariance matrix 
7, L is a lower triangle matrix. Given the nature of CD, we can parallel the decomposition. 

DGP is a multivariate Gaussian random vector. For a N dimension multivariate Gaussian 
random vector— i^^, ^, its covariance matrix 7 is defined as in terms of time correlation 
functions (kernel of influence functional), 

7 = (eT). (19) 

In order to generate the random vector, ^ in Eq. lAll with Cholesky decomposition, we 
need to generate A^ independent normal distributed random variables, Q, with Monte Carlo 
according to the normal distribution function, 2^exp(— -y)- The random vector ^ can be 
defined as, ^ = L * (^, where C = [Co; Ci) ' ' ' > Ca^]'^- With this construction, we can recover 
the following equality, 

(eD = {LVei^) = (LIL^) = 7, (20) 

where Cs = ^ since (OCj) = ^ij- Fo^' the complex- valued Gaussian process, since the 
covariance matrix of PDF in Eq. IB5I is not complex symmetric or Hermitian, finding the 
reliable algorithm of decomposition is challenging but possible^. We will discuss how to 
sample a complex-valued Gaussian process in the future paper. So previous work— 1^ shows 
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the convergence of the complex-valued Gaussian process is slow based their construction. 
But we anticipate that this construction with tweaks can improve the efficiency of the Monte 
Carlo for the complex-valued Gaussian process. In this paper, we will limit our method to 
the high temperature limit and real-valued Gaussian processed as used in Eq. [T3J 

B. Conditional Propagation Scheme 

Once we can sample the environment fluctuations $,{t) corresponding to influence func- 
tional, we can propagate the reduced density matrix according to Eq. [TH We propose a 
conditional propagation scheme to calculate the reduced density evolution by averaging over 
the Gaussian trajectories. In the propagation scheme, we need to generate one realization 
of the discrete Gaussian process ^ ffist using the CD method for the whole discrete time 
grid, [0,ti,t2, . ■ . , tAr„i, t] as discussed in the previous subsection. We can solve the stochas- 
tic Liouville equation in Eq. [13] by propagating the conditional density matrix using the 
fourth order Runge-Kutta scheme. We also can propagate the conditional density matrix 
with the following iterative scheme based on the time-sliced stochastic evolution operator 
exp[-tdt{Hs + Vm)], 

ps{t + dt\i{t + dt)) = exp [-2 dt {H, + V e(t))] p{t\i{t)) exp [i dt {H^ + V i\t))] . (21) 
The steps of the scheme can be described as: 

1. Generate ^ using the Monte Carlo according to PDF in Eq. IB5I using the CD method; 

2. Propagate the conditional density matrix ps{t\^(t)) to ps(t + dt\^(t + dt)) using Eq. [^ 
step by step. 

From the scheme, we can see that at every step, the density matrix is conditional on the 
environmental Gaussian random field, C,{t) and ^'{t). For the initial density matrix, ps{0) is 
conditional on the ^(0) so that we can write it as p5(0|^(0)). We choose the left end point, t, 
of the interval [t, t + dt] to define our stepwise evolution operator, exp {—i dt {Hs + V ^{t))) 
and exp [i dt {Hs + VC,'{t))). It might be interesting to explore the numerical scheme using 
the middle point or higher order approximation in the stepwise evolution operator to have 
better efficiency and accuracy at larger time step size, dt. 
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V. RESULTS AND DISCUSSION 

We introduce the conditional propagation scheme to compute the evolution of conditional 
reduced density matrix, ps(t\^) dependent on Gaussian trajectories, C,(t). Since the Gauss- 
Markov model has been used widely before, it can sever as a good benchmark model to 
validate the propagation scheme. Besides it, we will apply the scheme on the FMO system 
to study the dynamic effect of spectral densities. 

A. Benchmark with Gauss-Markov Model 

The Gauss-Markov Model has been extensively used to study exciton transport^. We 
take a symmetric dimer again as an example, which has the following stochastic Hamiltonian, 

2 

H = J2e,\k){k\ + J(|l)(2| + |2)(1|) + Se^\l){l\ + Se2\2){2\, (22) 

fc=i 

where the energy fluctuations, 5ei, are the Gauss-Markov process with exponential kernel 

{6e,{t)6e,{t')) = A^ exp{-^\t - t'\)6,,. (23) 

For the symmetric dimer (ei = 62), the local master equation (partial ordering prescription) 
of reduced density matrix is a set of coupled integro-differential equation^, 

^'^^^ - Uy) (24) 



dy 

-(1 - 2A^ g^{y)) p,{y) - 2A^ g^ ^'.(y), 



i's[y) _ /. oa2„/„\n„/„a oa2 



dy 

where P{t) = pu - P22{t), m = ^{p2i{t) - puit)), Ps{y) = P(y/2J), Mv) = Hyf^J), 

y = 2Jt, t is the time, and A^ = 27- Using the symmetric dimer, we can benchmark 
our conditional propagation scheme in two ways: 1. the first benchmarking, comparing 
the results of the conditional propagation scheme with the results of Eq. [2l] at the weak 
damping limit. A/ J ^ 1; 2. the second benchmarking, comparing the CD method with the 
Markov chain Monte Carlo (MCMC) sampling method which only works for the exponential 
kernel^. 

In the calculations, we choose ei = €2 = 0, J = 0.5, A = 0.05, and 7 = 1.0. Figure [H 
shows the first benchmarking results for the population pu — P22 and coherence P21 ~ Pu- 
The two results agree with each other very well. Figure [2] shows the second benchmarking 
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FIG. 1: Comparison of the results of conditional propagation scheme and the POP 
second-order master equation at the weak damping limit 

results for the population pn — P22 and coherence P21 — Pu- In the second benchmarking, 
the MCMC method is used to generate the independent Gauss-Markov processes C,{t). We 
have a good agreement as well. For both benchmarking, we use 5000 trajectories. 

B. Effect of Spectral Density in the FMO System 

Spectral density plays an important role in the dynamics of reduced system since the 
distribution of coupling coefficient decides how the energy flows into and out of the system. 
Mathematically, the relaxation process due to stochastic Gaussian environment is governed 
by the kernel of Gaussian process. The system of FMO has been studiecl^'^ extensively, we 
take the FMO model as an example to look at the influence of different kernels. The system 
Hamiltonian in the references^"— is used. The recent work- shows that the coherence in 
FMO system could depend on the initial condition. However, we in this paper will focus on 
how dephasing rates get manipulated by the spectral density given a fixed initial condition 
at site 1 and thereof coherence. We treat the Hamiltonian more as another model system. 

The real part of kernel 7(r — cr) is determined by spectral densities. 



LAr-a) 



dujJ{(jj) coth(/3/la;/2) cos(ci;t). 



(25) 



Ohmic spectral density has the form J{uj) = r]U) and it leads to the 6{t) spectral density^. 
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FIG. 2: Comparison of the results of the MCMC and CP methods (the imaginary part of 
the populations pn and P22 in Panels (a) and (b) are at the magnitude of l.Oe-10) 

However, Ohmic spectral Density with Lorentzian or exponential cutoff is used extensively 
in the study of non-Markovian excitation energy transfer in the light harvesting complex. 
Ohmic spectral density with Lorentzian cutoff or Drude spectral density is particularly 
popular because of its connection to Gauss-Markov process (exponential kernel) as revealed 
in Subsection IIIICI For the current model, the geometry of spectral density should play an 
role in the dynamics of reduced systems. For example, at the high temperature or strong 
coupling limit, the probability of multiphonon gets high and the phonon side band becomes 
closer to a Gaussian distribution^"^. The shift will change the reduced dynamics. The 
Gaussian spectral density is used to model the optical phonon band^S. It is interesting to 
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study the effect of a optical Gaussian phonon on EET— . It is also important to notice that 
the different cutoffs of the Ohmic spectral density will generate different non-6 kernels at 
hight temperature limit and different non-Markovian effect on the reduced system. 

In this section, we choose four different kernels: 1. the kernel of Drude spectral density, 
J{uj) = f]u ^, 2?i^2i in reference^; 2. the kernel of Ohmic spectral density with exponential 
cut-off J{ijj) = -qu^ exp{—u/uc) with Uc = 5ps~^; 3. the kernel of the first Gaussian 
spectral density rju \— exp(— ^'^"^7 ) with ujop = with cr = 7; 4. the kernel of the second 
Gaussian spectral density with a = 27. All the parameters are set up according to the 



paper-^. 



A2 



The Drude spectral density is used as the reference reorganization energy 77 = —^. In 
order to compare results at the same level, we normalize the four spectral densities to the 
reference reorganization energy, J^ dujJ{u)/uM: The superimposing of the four different 
spectral densities is shown in Figure |3l In the figure, ohmic with exponential cut has the peak 
around 15cm~^, the Drude and first Gaussian spectral density around 25cm~^, the second 
Gaussian spectral density around 40 cm~^. The major difference between the Drude and 
first Gaussian is the weight of high frequencies. The Drude has longer tail and more weight 
at the high frequency domain. Ohmic spectral density with exponential cutoff, compared to 
other spectral density, mostly concentrates in the low frequency domain. The comparisons 
of the four corresponding kernels are shown in Fig. HJ Panel (a) in Fig. H] shows that the real 
part of the kernel of the Drude spectral density is almost identical to the exponential kernel 
and tells us that at high temperature, the model with Drude Ohmic spectral density will 
collapse to the Gauss-Markov Model. We also can observe that even though the Drude and 
first Gaussian spectral densities have the same peak, they have different kernels in Fig. HI 
The real part of the kernel of the Drude spectral density has longer curve and bigger area 
under the curve compared to the first Gaussian spectral density. 

Fig [5] shows population dynamics pii{t) for the four different spectral densities. Fig [6] 
clearly shows the imaginary parts of populations converge to zero. We can see that the Drude 
spectral density and second Gaussian spectral density give the fastest decay since both have 
more weight in the high frequency domain and smaller area under the real part of Kernel 
curve in Figure HI Our conditional propagation scheme uses 5000 sampling trajectories for 
all the FMO results. 

We also can see that while the Drude and first Gaussian spectral densities have the same 

15 



peak, they have different areas under kernel curves and different decay rate. Based on 
these observations, we can draw a simple conclusion that the low frequencies lead to the 
non-Markovian effect (memory effect) and high frequencies are more associated damping 
and decayi^. But the areas under the kernel curves essentially decide decay/damping rates. 
However, the connection between spectral density and the area under the kernel curve is non- 
linear. It suggests that tweaking with the shape of spectral density give us more interesting 
insight to the interplay between coherent and incoherent motions. 



Spectral Density 



Drude Ohmic 
Exp-Cut 
Gaussian 1 
Gaussian 2 




100 150 

Freq cm^^ 



250 



FIG. 3: We present four spectral densities: Ohmic spectral density with Lorentzian cutoff, 
Ohmic spectral density with exponential cutoff, and spectral densities with Gaussian cut 

off with cr = 7 and a = 27 (7 = 5 ps~^). 

Figure [7] shows that the comparisons of the population dynamics for three different Ohmic 
spectral densities with three different exponential cutoffs, Uc = 2,5,10 ps~^. It is clearly 
shown that by fixing the reorganization energy, the more the spectral density is shifted to 
the higher frequency, the faster the relaxation/decay is. 

In summary, based on the model we have here, the low frequencies is associated the 
memory (non-Markovian) effect and high frequencies are more associated with damping or 
decay. When we change the weight among low and high frequencies, we can change the 
curve of kernel and the relaxation process. All these arguments are based on the fixed initial 
conditions. 
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FIG. 4: Kernels of four different spectral densities, Lorentzian (a), Exponential Cut (b), 
Gaussian with a = 7 (c) and Gaussian with a = 27 (d) (7 = 5 ps~^) 

VI. CONCLUDING REMARK 

In the paper, we review coherent state path integral and influence functional. By exploit- 
ing the similarity between GCF and influence functional, we can construct complex-valued 
Gaussian processes to deconvolute the dynamics of reduced system temporarily. In the con- 
struction, we build the covariance matrix of PDF of discrete Gaussian propose according to 
convolution kernel of influence functional. Using the CD method, we can sample the dis- 
crete Gaussian process. On top of it, we propose conditional propagation scheme to simulate 
dynamics of the reduced system under the influence of different spectral densities. 

Using the CD method and conditional propagation scheme, we examine EET in the FMO 
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FIG. 5: Population Dynamics of FMO, pui^t), in the site representation for the kernels of 
four different spectral densities: Lorentzian, Exponential Cutoff, and Gaussians with two 

different a's. (a) shows the results for Lorentzian, (b) the results for the kernel of 
exponential cut, (c) and (d) the results of kernels of Guassians with two different widths, 

(7 = 7 and cr = 27 (7 = 5 ps~^) 

complex under the influence of different spectral densities. We find that the geometry of 
spectral density can change the reduced system dynamics. In this paper, we find that the 
the low frequencies is more associated with the non-Markovian effect (memory effect) and 
high frequencies are more associated damping and decay. Since the connection between 
spectral density and the area under the kernel curve is nonlinear, tweaking with the shape 
of spectral density may give us more interesting insight to the interplay between coherent 
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FIG. 6: Corresponding imaginary part (in reference to Fig [5]) of the populations, Pii{t) 
(y-axis scale is at the magnitude of l.Oe-9 or l.Oe-10 for panel (d)) 

and incoherent motions. 

For the future work, we will extend the CD method to sample complex-valued Gaussian 
process according to the complex-valued non-Hermitian covariance matrix. We also like to 
apply the method to study of the correlated initial condition for the open quantum system. 
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FIG. 7: Panel (a) shows the spectral densities with differential exponential cutoffs, 
Wc = 2, 5, 10 ps~^ and Panels (b), (c), and (d) show the comparison of the results, the 
population dynamics of reduced systems, pu, for the three different spectral densities. 
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Appendix A: DGP and GCF 

This representation of the stochastic DGP, ^{t) is characterized by its kernel, i.e. two 
time correlation function, 7(t) = (^(t)^(O)). In this section, we would like to show how 
we can construct the proper Gaussian process by treating the influence functional as a 
GCF. Physicists are more familiar with the Fokker-Plank equation when we discuss the 
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Brownian motion. From the practical perspective, Gaussian fluctuation's dynamic properties 
are determined by the correlation. Most of dynamical systems that we meet are stationary, 
which means the two-time correlation function is only dependent on the time difference, 

DGP can be deflned as a random vector ^ on the discrete time lattice, [to, ■ ■ -tn], defined 
as, 

i=mo),m)r--,atn-l),atn)r, (Al) 

where to = and t„ = t. For the simplification of notations, we replace $,{ti) with ^i. 
The joint probability density function (PDF) of DGP ,^ is a multivariate Gaussian function. 
To define the multivariate Gaussian PDF, we need two structure parameters, mean vector, 
m = [mi, 7712,- ■ ■ ,?Tin]^ where rrii = {C,i), and covariance matrix, 7 = [-fij] where 'jij = 
{{C,i—mi){^j—mj)). Based on the mathematical properties of multivariate Gaussian function, 
PDF of DGP is defined as, 

P(0 = Nexp[-{i - myj-\i - m)], (A2) 



where 



7 



((^0 - "io)(^o - "^0)) ((^0 - "io)(6 - "^1)) 
((6 - ?Tii)(^o - "^o)) ((6 - "^i)(6 - "^1)) 



(A3) 



and A^ = , .„L 11/2 the normalization factor and | ■ | determinant. Without loss of generality, 
we assume that m = in the context of mean field. 

The Gaussian process C,{t) is the continuous limit to the DGP, C, when dt — )■ 0. What is the 
meaning of the limitation? How do we express the limitation in terms of what quantity. PDF 
of DGP ^ governs the discrete Gaussian environmental fiuctuations. However, PDF's for 
continuous Gaussian processes don't exist. In other words, we can't take the continuous limit 
on PDF. The Fourier transform of PDF's, i.e. GCF, provides the equivalent information to 
PDF's. We can use GCF as the quantity to define the continuous limitation of DGP. In the 
following subsection, we will define the GCF and show the connection of GCF to infiuence 
functional. 

Corresponding to Eq. IA2t the GCF of DGO is defined as, 

(exp(-z Y, ^^^" dt))^ = exp(-^ ^ F7S dt^), (A4) 



2 
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where we assume homogeneous time lattice, ti — tj_i = dt, and s = 
[s(to),s{ti),s{t2), ■ ■ ■ ,s{tn)]\ the dummy vector (a discrete deterministic process). This 
is essentially the average of the Fourier Transform of PDF. By taking t/t — )■ in Eq. IA4t 
the integral sign, J, will replace the sum sign, ^. Given the symmetry of the kernel, 
7(— t) = 7(t), we have the continuous-time Gaussian process GCF, 

(exp(-i Jl^s^idt)) = exp[-i/grfr/Qdas(r)7(r,a)s((T)] 
= exp[-/Q rfr/([c/(T s(r)7(r,(T)s(a)]. 
The continuous Gaussian process kernel, 7(0", r) is the the continuous limit of the discrete 
covariance matrix 7 in Eq. IA3I 

GCF for real (continuous-time) Gaussian process will give us the lead to construct the 
Quantum Gaussian random field to reproduce the influence functional. The simple compar- 
ison already shows that influence functional is similar to GCF. Therefore, constructing the 
proper covariance matrix, 7, for the complex-valued Quantum Gaussian random field, we 
can achieve our goal to map the convolution of influence functional to a random field. We 
will show how to construct the mapping in the next appendix section. 

Appendix B: Influence Functional and Complex- Value Gaussian Process 

In order to draw the linkage between influence functional and general characteristic func- 
tion, we need to discretize influence function. Given that influence function in Eq. |H] is a 
time ordered double integral and the symmetry of the Kernel 7(— t) = 7^(t), we have 



F{Q,Q';tf,U) = exp{- dr da {V{Q{r)) - V{Q'{t)) (Bl) 

Jti Jti 

[7(r - a)V{Q{a)) - ^^{r - a)V{Q\a))]} 

= expif^dr f da -V{Q{r)y/{T-a)V{Q{a)} 
Jti Jti 

+ expif'dr f da - V (Q' (t))^^ {r - a)V {Q' (a)} 
Jti Jti 

+ exp{f'dr f da +V{Q'{T))^ir-a)V{Qia)} 
Jti Jti 

+ expifdr f da +V{Q{T))^\T-a)V{Q'{a)} 
Jti Jti 

We define N homogenous discrete time grid tj where i = 0,1,2,- ■■ ,N, to = ti, tw = tf, 
and dt = {tj — ti)/N , Vt = V{Q{ti))dt, VI = —V{Q'{ti))dt and covariance matrix element, 
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lij ~ 7(^« ^ ^3)1 Tij ~ iH'^i ~ '^j)- The discrete version of F{Q, Q'; tf, ti) can be defined as, 



It] 

N i 



N i N i N i \ 

-p ( - E E ^^^.^ - E E ^'^ii^ - E E ^ 7.. V, - E E ^^l^ (^2) 

\ j=o j=0 i=0 j=0 i=0 j=0 i=0 i=0 / 

Following the constructive approach in Sec. Illlt the covariance matrix for discrete 
complex- value Gaussian process, ^(t), is defined as. 



7 



(ei^o) (Cl^l) 






(B3) 



where 



(B4) 



Sn 

where ^j = ,^(tj) and ,^j' = ^'(ti) and the means of ^j and ^j' are zeros, {^iC,j) = {C,jC,i) = 'Jij, 
(C'e;) = (e^:) = ll W,) = (e;e.) = ll where ^ > j and Uj) = (e^.) = 7^j where 
2 < J. Also 7ii = 7!^. Therefore {^S) = {^-Q = (d^i) = (^j^-) = 7^ = 7!^. This matrix is 
non-Hermitian which is one intrinsic property of quantum open systems. 
The corresponding complex-valued PDF of ^ can be expressed as. 



where A^ is the normalization factor. The corresponding GCF is defined as, 

(exp(-z ^V^^))-=exp(-^V^7V), 



(B5) 



(B6) 
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where the dummy vector, 



V 



Vo 

Vn 



v 

After comparing Eq. IB6I and Eq. IB21 we can found our construction of complex- value Gaus- 
sian noises generates GCF equal to influence functional for the discrete version. In order to 
prove the equality, we define the follow two vectors Vi and V2 and four blocks 'Jn, 722) 
7i2, and 731: 



Vi 



Vi 



Vi 



N 



111 



7l2 



{i^Q (66) 

(66) (66) 



^0 
Vi 



V' 



722 



,721 



(6^0) (6^i) 
(66) (ei6) 

(66) (66) 
(66) (66) 



Eq. 



can be re-written as, 



exp(-lv^7iiVi - \vl-t^,V2 - \vll,,Vi - ^^^7121^2) 



(B7) 



(B8) 



(B9) 
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" T 

Since 'jn and 722 ^^^ symmetric matrices according to their definitions, exp(— ^ V^ 7ii^i) = 
exp(- E!Io Ej=o '^i7ijV,) and exp(-|V2 722^^2) = exp(- J2f=o E}=o '^/^I,'^/)- For the re- 
maining two terms, we have to consider the following eqaulity, 

exp ( - ^^^712^^2 - ^V^72iV^i) = (BIO) 

N i , N N 



i=0 j=0 i=0 j=i+l 

. N j N N 

-5 E E i^M v/ - 5 E E '^M'^/) = 

j=o j=o j=o «=i+i 

^ N i N N 

^^p ( - 9 E E ^1 ^i^^2^' - 9 E E ^i^i^^2 



9 / ^ / J • L I Li • i cy 

i=0 j=0 i=0 J=i+1 

4EEt^'^at'i-^EE''.'^S/,K,o = 

j=0 i=0 j=0 i=i+l 

^ AT J AT AT 

^^p ( - 2 E E ^i^i^i - 2 E E ^i^i^^2 

i=o i=o j=o i=j+i 

^ AT 4 Af AT 

-sEE^-^^i^-sEE^^^ii^)- 

i=0 j=0 j=0 i=i+l 

If 7*^2 = 721 = lij when 2 > j and 7^-'2 = 721 = 7ij when i < j which are satisfied in the 
covariance matrix Eq. IB3t the equality becomes, 

^ N i ^ N N 

^^p ( - 2 E E ^^^^ - 2 E E ^i^i^^2 (Bii) 

1=0 j=0 i=0 j=i+l 

N i N N 

-2 E E ^^^i^ - 2 E E ^^^1^) = 

i=0 j=0 j=0 j=i+l 

N i N N 

^^p ( - E E ^i^^i^^2^' - E E ^^^) ■ 

i=0 j=0 i=0 j=i+l 

It is easily to show that X]i=o Si=i+i ^2'72i^i = J2i-o J2] ^0^2 121^1 ■ With this, we com- 
plete our proof, 

(expH^V^O)^= (B12) 

i 
/ N i N i N i N i N 

-p - E E f-^u^. - E E ''i-i^i - E E ^^f. - E E "-in' 



4=0 j=0 i=0 i=0 i=0 jr=0 i=0 j=0 
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because V-l corresponds Vi and V2 corresponds to V-. By taking the continuous limit of 
Eq. IB6t the corresponding continuous-time GCF should match the influence functional in 
Eq. [HI In other words, we prove the equality in Eq. O 



Appendix C: Complex Gaussian and Classical Gaussian Noises 



At the classical high temperature, the complex covariance matrix will become 



7 



(aeo) (aei) 






(Cl) 



which becomes a symmetric real-valued matrix, i.e., the four matrices become equal to each 
other, 7]^^ = 722 = 721 = 7i2- Iii order to reproduce the high temperature real- valued 
covariance matrix, the random vector ^ can be simplified to be, 

6 



Co 

6 



(C2) 



z.e. e: = e.. 
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Appendix D: Gaussian Process for Dimers 

The influence functional can be defined as (for example in the reierenc^^] 



I{s^, st,...,s,SQ,s^,...,s') = (Dl) 

rYj.l^-i{Hi{s,X)+IxHB{X))dt/2^-i{Hi{s+_-^,X)+IxHB{X))dt 
^-i{Hi{s+ ,X)+IxHB{X))dt^-i(Hi{s+ ,X)+IxHB(X))dt p,^-. 
^iiHi(s- ,X)+IxHB(X))dt^i{Hi{s- ,X)+IxHB{X))dt 
^i{Hiis^_^+IxHB{X)),X)dt^iiHi{s',X)+IxHB{X))dt/2i 

where Si is the discrete system variables, since 

f Hn \ 
Hi = VxX = \ (D2) 

\ Hj,J 

is diagonal, where Hji = Vi * X and if/2 = V2 in the dimer Hamiltonian in Eq. [151 

^-^iHJis,,X)+IxHB(X))dt ^ f exp(-2(yi(s,) * X + HB{X))dt) 

y eM-i{V2{si)*X + HB{X))dt). 

(D3) 

And the following influence functional, 

7(4, s]^, • • • , s. So , sj", . . . , s') = (D4) 

^-iHei{s,X)dt^-iH^i{s+^^,X)dt ^ ^ ^ ^-iHei{s+ ,X)dt q 

Q g-iHe:2(s,X)dt^-iHe2[sti_-t_,X)dt ^-iHe2{s^ ,X)dt 

xpb(O) 

^.H ei{s- ,X)dt _ _ _ ^iH^-i(s-^_^,X)dt ^iHel(s' ,X)dt Q 

Q ^iHe2{s~ ,X)dt ^iHe2{s]^_j,X)dt^iHe2{s',X)dt 

where pb(0) = exp(-/3i7,,)/Tr(exp(-/3ifB)), H,, = V^{s,) * X + //^(A) and iJ,2 = V^2(s«) * 
A + Hb{X). After integrating over the degree freedom of bath for the diagonal matrix 
elements in Eq. ID41 we get two separate influence functionals for each diagonal matrix 
element. As a result, we should have two independent Gaussian random processes for each 



Tr 
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